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Final  Report 

Grant  #  W911NF-09-1-0199 

Ab  Initio- Based  Predictions  of  Hydrocarbon  Combustion  Chemistry 

PI:  Donald  L.  Thompson 

Department  of  Chemistry,  University  of  Missouri-Columbia,  Columbia,  MO  65211 

Overview.  The  practical  objective  of  the  research  was  to  address  some  of  the  challenges  to 
developing  improved  hydrocarbon  combustion  models.  This  requires  new  theoretical  and 
computational  methods  for  predicting  reaction  rates,  particularly  for  extreme  pressures.  We  have 
developed  and  applied  a  novel  approach  for  computing  the  effects  of  pressure  on  reactions  and 
relaxation  of  radicals  and  molecules,  which  we  have  used  in  studies  of  the  relaxation  of  excited 
species  in  an  Ar  bath.  This  method  allows  us  to  avoid  (and  evaluate)  the  usual  assumption  of 
bimolecular  collisions,  which  breaks  down  at  high  pressures.  For  predictability  we  need  to 
develop  efficient  methods  for  using  ab  initio  potential  energy  surfaces  (PESs)  computed  with 
high  levels  of  quantum  chemistry  theory  to  provide  accurate  rates  of  elementary  reactions.  The 
methods  all  make  use  of  classical  molecular  dynamics. 

The  Influence  of  Pressure  on  the  Relaxation  of  Excited  Radicals  and  Molecules.  A  major 
thrust  of  the  research  was  on  the  effects  of  pressure  on  the  relaxation  of  excited  species  in  the  gas 
phase.  Our  main  interest  was  to  explore  ranges  of  pressure  where  the  relaxation  occurs  primarily 
in  “isolated”  bimolecular  collisions  to  levels  where  the  bimolecular  in  no  longer  valid.  The 
relaxation  of  radicals  and  molecules  is  usually  assumed  to  occur  in  collisions  of  the  excited 
species  with  a  single  other  species.  Clearly  this  bimolecular  mechanism  must  cease  to  be  valid  at 
extremely  high  pressures.  Our  interest  was  to  determine  the  change  from  the  bimolecular  process 
to  relaxation  processes  in  which  more  than  one  other  molecule  or  radical  in  involved. 

Molecular  dynamics  simulation  of  the  relaxation  of  three  different  species  (CH3NO2,  HO2, 
and  OH)  in  an  argon  bath  gas  at  pressures  over  the  range  10  atm  to  400  atm  and  at  temperatures 
of  300  K  (CH3NO2  and  OH)  or  800  K  (HO2)  were  carried  out.  The  three  species  were  chosen  for 
their  considerable  differences  in  size,  and  molecular  and  vibrational  structures.  Acceptably 
accurate  intramolecular  potential  energy  surfaces  were  available  for  all  three.  The  interactions  of 
CH3NO2  and  HO2  with  Ar  were  approximated  by  pairwise  additive  Buckingham  (exp-6) 
potentials  with  the  atoms  in  the  excited  species  treated  as  the  nearest  rage  gas,  and  the 
parameters  determined  by  using  combination  rules  and  parameter  values  available  in  the 
literature;  which  has  been  shown  to  provide  reasonably  accurate  intermolecular  forces  between 
molecules  and  radicals.  The  Ar-H  interaction  in  OH  was  represented  by  the  Lennard-Jones 
potential  and  the  Ar-O  as  in  the  other  species,  that  is,  with  the  exp-6. 

Classical  molecular  dynamics  simulations  were  carried  out  for  1000  Ar  atoms  and  one 
initially  excited  species  in  a  cell  with  periodic  boundary  conditions.  The  molecule  or  radical  was 
instantaneously  excited  by  either  statistically  distributing  energy  among  the  internal  degrees  of 
freedom  (50  kcal/mol  for  CH3NO2  and  45  kcal/mol  for  HO2),  but  in  the  case  of  OH  the 
relaxation  of  the  specific  state  v=2,  J=0  was  studied.  The  translational  energies  of  the  excited 
species  and  Ar  atoms  were  selected  from  the  Maxwell-Boltzmann  distribution. 

A  comparison  of  the  relative  decay  rates  for  rotation  and  vibration  at  substantial  pressure  is 
shown  in  Fig.  1  for  the  three  systems.  In  all  three  cases,  as  expected,  rotational  relaxation  is  one 
to  ten  orders  of  magnitude  faster  than  vibrational  relaxation.  The  results  clearly  show  the 


1 


equilibration  of  the  rotational  energy  to  values  approaching  1.5 kT  for  CH3NO2  and  HO2  (~0.9 
kcal/mol  (T  =  300  K)  and  -2.4  kcal/mol  (T  =  800  K),  respectively,  and  kT  for  OH  (-0.6 
kcal/mol).  However,  the  vibrational  relaxation  rates  vary  among  the  three  systems  by  at  least  an 
order  of  magnitude;  the  rate  for  CH3NO2  is  the  fastest  and  for  OH  by  far  the  slowest.  This 
variation  in  relaxation  is  likely  due  to  the  different  types  of  vibrational  motion  possible  in  these 
species;  for  example,  CH3NO2  has  a  low  frequency  hindered  rotor  that  can  couple  well  with  the 
motions  of  the  bath  gas  atoms. 
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Figure  1:  The  ensemble  averaged  vibrational  energy  (in  red),  rotational  energy  (in  blue),  and 
translational  energy  (in  green)  as  a  function  of  time  for  CH3NO2  in  a  (300K,  100  atm)  Ar  bath 
(left),  for  HO2  in  a  (800K,  100  atm)  Ar  bath  (middle),  and  OH  in  a  (300K,  50  atm)  Ar  bath 
(right). 


Figure  2  shows  the  normalized  excess  vibrational  energy  Enrm-vib(t)  for  all  the  simulations. 
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Figure  2.  The  normalized  excess  vibrational  energy  (see  text)  as  a  function  of  time  for  CH3NO2 
in  a  300K  Ar  bath  gas  at  eight  pressures  (left),  for  HO2  in  a  800K  Ar  bath  gas  at  five  pressures 
(middle),  and  for  OH  (v=2,J=0)  in  a  300K  Ar  bath  gas  at  one  pressure.  The  solid  lines  are  the 
simulation  results  while  the  broken  lines  are  the  LS  fit  (see  text). 
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Enrm-vib(t),  the  vibrational  energy  above  the  thermal  vibrational  energy,  is  plotted  as  a  function  of 
time,  including  the  temperature  increase  of  the  bath  gas  due  to  molecular  energy  lost  from  the 
initially  excited  molecule  up  to  that  time.  The  results  are  normalized  to  the  excess  energy 
appropriate  to  the  initial  absolute  value  of  the  energy.  The  normalization  allows  the  different 
systems  and  pressures  to  be  conveniently  compared.  The  dashed  curves  in  Fig.  2  are  fits  of  the 
trajectory  results  to  the  Lendvay  and  Schatz  (LS)  [G.  Lendvay  and  G.  C.  Schatz,  J.  Phys.  Chem. 
95,  8748  (1991).]  equation: 

This  functional  form  for  the  decay  of  the  normalized  energy  Emm(t )  has  the  virtue  that  it 
involves  only  two  parameters,  and  each  parameter  responds  to  different  temporal  regions  of  the 
decay.  At  /=0,  the  ki  parameter  equals  dEnnn/dt,  the  initial  rate  of  decay.  To  the  extent  that  this 
form  is  appropriate,  k,  will  be  determined  primarily  by  the  early  behavior  of  the  decay  curve.  The 
parameter  m  determines  curvature.  When  m  =  1,  Eq.  (1)  can  be  shown  to  collapse  to  a  single 
exponential  function  with  rate  ki.  For  m  >  1,  Emm(t)  rises  up  from  the  single-exponential  straight- 
line  behavior  on  a  semi-log  plot,  indicating  that  as  excitation  level  of  the  species  decreases  the 
relaxation  rate  decreases.  For  m  <  1  (including  negative  values),  Emm(t)  curves  down  from 
single-exponential  (linear)  behavior,  indicating  that  as  excitation  level  decreases  the  relaxation 
rate  increases.  Figure  2  shows  positive  curvature  at  all  pressures.  For  the  lowest  pressure,  10  atm, 
it  appears  that  a  single  decay  rate  applies  over  1000  ps  for  CH3NO2  and  HO2;  however,  as  the 
results  in  the  inset  show,  when  the  process  is  followed  for  5000  ps,  there  is  clearly  positive 
curvature.  While  the  LS  functional  form  fits  the  results  quite  well,  we  see  it  as  merely  an 
empirical  fit  without  a  theoretical  basis  or  justification.  We  were  unsuccessful  in  finding  a 
physical  explanation  for  bi-exponential  behavior,  something  that  has  long  been  observed,  used  in 
modeling  of  combustion  and  other  complex  processes,  yet  never  shown  to  be  more  than  an 
empirical  fitting  function. 


Figure  3.  The  optimized  values  of  m  (in  red  and  corresponding  to  the  red  axis)  and  ki  (in  blue  and 
corresponding  to  the  blue  axis)  as  functions  of  density.  The  values  are  from  fits  to  the  LS 
equation,  shown  in  Fig.  2.  The  densities  correspond  to  the  pressures  in  Fig.  2. 
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Plots  of  the  parameter  values  m  and  as  functions  of  density  are  shown  in  Fig.  3.  The  density  is 
determined  by  the  size  of  the  periodic  box  in  the  simulations  at  each  pressure  and  temperature 
(this  accounts  for  any  deviation  at  higher  pressures  from  the  ideal  gas  law).  The  different  density 
scale  for  HO2  reflects  the  fact  that  at  nearly  three  times  the  temperature,  density  over  the  same 
pressure  range  is  about  1/3  less. 

As  one  can  see  from  the  red  points  and  curves  in  Fig.  3  there  is  positive  curvature  in  all  cases 
with  the  HO2  curve  considerably  more  curved  than  that  for  the  other  two  systems.  The  solid 
curve  in  Fig.  3  is  a  fit  of  the  density  dependence  of  m  to  the  functional  form  a  +  bpecP ,  where  p 
is  the  density.  This  fit  shows  that  while  both  CH3NO2  and  HO2  show  some  pressure  dependence 
to  m,  the  m  values  lie  in  a  relatively  narrow  band.  The  persistence  of  this  band  of  curvature  over 
such  a  wide  variation  of  conditions  suggests  that  the  internal  structure  of  molecule  is  involved  in 
the  relaxation  process.  For  OH,  the  simplicity  of  this  internal  structure  permits  a  deeper 
interpretation  of  the  origin  of  curvature.  OH  has  only  one  vibrational  degree  of  freedom  and  Fig. 
1  shows  that  only  minor  amounts  of  vibrational-rotational  coupling  occur  for  the  initial  OH  state 
of  (v=2,J=0).  One  can  then  model  the  vibrational  energy  decay  in  terms  of  decay  rates  from 
specific  vibrational  levels  and  obtain  an  analytic  solution.  Because  the  energy  gap  between  OH 
vibrational  levels  is  so  large  (3737  cm'1  fundamental  frequency),  the  only  substantial  rates  are 
downward  rates  to  adjacent  levels,  i.e.,  kv,v.  1  where  v  is  the  vibrational  quantum  number  of  the 
initial  state.  From  perturbation  theory  for  harmonic  oscillators,  it  is  well  know  that  kv,v-i  is 
proportional  to  v.  In  other  words,  as  energy 
leaks  out  of  OH  in  to  the  bath  gas,  the 
increasing  populations  in  lower  OH 

vibrational  levels  experience  slower  rates  of 
further  de-excitation.  Nonetheless, 

substitution  of  kv<v.  1  =  vkio  into  the  analytic 
solution  produces  a  single  exponential  decay 
(m  =  1)  given  by  k\o.  The  only  way  to  get  m 
>  1  is  for  kVtV- 1  >  vk\o.  This  can  be  achieved 
by  recognizing  that  anharmonicity  of  the 
vibrational  levels  causes  the  energy  gap 
between  adjacent  levels  to  decrease  as  the 
energy  of  the  oscillator  increases.  Kohno  et 
al.  [N.  Kohno,  J.  Yamashita,  C.  Kadochiku, 

H.  Kohguchi,  and  K.  Yamasaki,  J.  Phys. 

Chem.  A  117,  3253  (2013).]  in  an  OH/He 
single  collision  experimental  study  observed 
that  kv>v. i/v  depended  exponentially  on  the 
anharmonic  energy  gap.  Inserting  an 
exponential-gap  dependence  in  the  analytic 
kinetics  model  produces  positive  curvature, 
as  is  found  in  the  simulations.  This  suggests 
that  the  positive  curvature  present  in  all  the 
results  in  Fig.  2  are  due  at  least  in  part  to  the 
reduced  effects  of  anharmonicity  at  lower 
internal  energies  within  the  three  molecules. 
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The  results  in  Fig.  3  show  that  the  ki  parameter  assumes  values  that  increase  monotonically 
over  the  factor  of  ~40  variation  in  the  density  for  both  CH3NO2  and  HO2.  However,  the  results  in 
the  figure  also  clearly  show  that  ki  is  not  directly  proportional  to  density  over  the  entire  interval. 
At  a  pressure  of  -100  atm  (density  -0.002  A'3  for  CH3NO2  and  -0.001  A'3  for  HO2),  ki  begins  to 
drop  below  what  is  expected  from  a  direct  proportionality.  Recall  that  ki  is  the  relaxation  rate  at 
t= 0  for  the  LS  fit.  The  initial  conditions  of  the  molecule  and  radical  are  statistically  identical  for 
all  densities.  This  implies  that  the  loss  of  direct  proportionality  is  a  consequence  of  a  change  in 
the  nature  of  the  interaction  of  Ar  bath  with  the  excited  species  as  a  function  of  density.  The 
solid  curve  in  Fig.  4  is  a  fit  of  the  density  dependence  of  ki  to  the  functional  form 

=  Pp(l)p/ti  +  TjL>  1  Pp(L)p(L/c2), 

where  p  is  the  density  and  Pp(L)  is  the  normalized  probability  of  having  L  atoms  colliding 
simultaneously  with  the  molecule,  and  kt  and  k2  are  adjustable  molecular  relaxation  rates  per  Ar 
colliding  atom.  A  simple  combinatorial  model  based  on  the  size  of  the  molecule  and  the  size  of 

the  Ar  atom  leads  to  an  analytic  expression  for  Pp(L)  and  suggests  that  at  the  highest  pressure 
(400  atm)  single  collisions  are  rare  and  simultaneous  collisions  involving  -5  Ar  atoms  is  the 
norm.  Roll-off  of  the  initial  rate  from  the  low-density  extrapolation  occurs  because  the  rate  for 
collision  events  with  L  Ar  atoms  grows  with  L  more  slowly  than  L  times  the  rate  for  collision 
events  with  one  Ar  atom. 

The  relaxation  of  normalized  excess  rotational  energy  can  be  analyzed  with  the  LS  fit  in  a 
manner  similar  to  Fig.  2  for  CH3NO2  and  HO2.  The  resulting  m  and  ki  parameters  are  plotted  in 
Fig.  5  in  a  manner  similar  to 
Fig.  3  for  vibrational  energy. 

There  are  two  significant 
differences  between  Fig.  5 
for  rotational  energy  and  Fig. 

3  for  vibrational  energy. 

First,  the  m  values  for 
rotation  are  relatively 
similar  for  the  two  species 
but  a  factor  of  two  different 
for  vibration.  Second,  for 
both  species  the  rotational  ki 
is  directly  proportional  to 
density  over  the  interval  of 
simulated  densities;  the 
deviation  from  direct 
proportionality  at  higher 
pressures  does  not  occur. 

The  reason  for  this  is  the 
subject  of  on-going  study. 

This  work  provides  further  insight  into  the  influence  of  increasing  pressure  on  relaxation  of 
internally  excited  species,  however,  it  is  clear  that  a  full  understanding  will  require  further  study. 
The  present  study  shows  the  breakdown  in  the  bimolecular  collision  assumption,  but  does  not 
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Figure  5.  The  optimized  value  of  m  (in  red  and  corresponding  to 
the  red  axis)  and  k,  (in  blue  and  corresponding  to  the  blue  axis)  as 
a  function  of  density  for  LS  fits  to  normalized  excess  rotational 
energy.  The  CH3NO2  results  (left)  and  the  HO2  results  (right) 
have  densities  corresponding  to  the  pressures  listed  in  Fig.  2. 
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provide  sufficient  knowledge  of  the  mechanisms  to  allow  us  to  explain  the  bi-exponential 
behavior  of  the  decay  rates.  Thus,  while  we  have  advanced  our  understanding  of  pressure  effects, 
further  work  is  needed  to  achieve  sufficient  knowledge  of  the  effects  to  either  devise  empirical 
relationships  with  wide  applicability  (e.g.,  to  cover  the  many  species  involved  in  combustion)  or 
to  have  a  satisfactory  theoretical  (including  mathematical)  model. 

Accurate  PESs.  Our  initial  focus  was  on  the  development  of  ab  initio- based  PESs  for  use 
in  MD  simulations  of  reactions.  We  began  by  considering  the  lower  cost  electronic 
structure  methods  for  C2H5.  We  first  investigated  empirical  potentials  especially  designed 
for  the  study  of  hydrocarbon  systems,  namely  Reaxff  and  AIREBO.  Both  are  cheap  to 
compute  and  are  close  to  acceptable  accuracy  for  C2H5.  We  also  examined  Density 
functional  theory  (DFT),  which  are  quickly  being  improved.  Some  give  results  comparable 
to  commonly  used  ab  initio  methods  at  considerably  lower  computational  cost.  We  tested 
the  commonly  used  functionals  B3LYP,  PBE,  and  PW91;  and  found  them  not  to  be 
sufficiently  accurate.  We  found  that  using  the  aug-cc-pVTZ  basis  set  that  the  M05 
functional  gives  acceptable  results.  We  are  currently  using  the  M06,  which  is  presumably 
more  accurate  than  M05,  to  develop  automatic  procedures  for  generating  global  PESs  and 
for  performing  direct  dynamics  simulations.  We  also  carried  out  direct  dynamics 
simulations  of  the  C2H5  dissociation  and  isomerization  using  DFT  M05-2X/CBSB7.  The 
isomerization  and  dissociation  reaction  barriers  are  almost  the  same  height,  thus  a  good 
test  case.  Direct  dynamics  simulations  for  a  system  of  this  size  with  the  high  accuracy  we 
require  are  expensive.  The  average  propagation  time  for  a  trajectory  at  150  kcal/mol  is 
about  two  weeks.  While  this  makes  the  calculations  feasible  (for  this  size  radical)  with  the 
use  of  clusters,  we  wish  to  lower  the  computational  cost.  We  will  continue  exploring  ways 
of  significantly  reducing  the  cpu  time  for  these  kinds  of  calculations. 

We  have  also  searched  for  suitable  reference  potentials  that  we  can  use  within  a 
strategy  we  developed  a  few  years  ago  in  which  we  fit  the  difference  between  an 
approximate  PES  and  ab  initio  data.  The  motivation  for  this  approach  is  that  the  difference 
can  be  considerably  "flatter”  than  the  actual  potential,  thus  requiring  many  fewer  ab  initio 
points  for  an  accurate  representation  of  the  ab  initio  PES.  our  focus  has  been  on  C2H5.  We 
have  investigated  ReaxFF  and  AIREBO.  Both  are  cheap  to  compute  and  are  close  to 
acceptable  accuracy  for  C2H5.  We  have  used  the  aug-cc-pVTZ  basis  set  with  the  M05 
functional  to  generate  10,000  data  points  for  checking  out  fitting  scheme.  If  an  effective 
practical  method  is  applicable  to  this  radical  it  could  be  used  for  other  larger  systems.  The 
isomerization  (the  transfer  of  a  H-atom  from  the  CH3  end  to  the  CH2  end)  and  dissociation 
reaction  barriers  are  almost  the  same  height.  These  competitive  reaction  channels 
represent  an  interesting  fundamental  problem  about  which  little  is  known,  thus  we  have, 
as  we  work  to  develop  the  IMLS-based  methods,  also  made  a  study  of  the  effects  of 
isomerization  on  dissociation.  Knowing  how  such  an  isomerization  affects  the  rate  of 
dissociation  could  provide  insights;  for  example,  different  energy  distributions  in  products 
that  might  affect  their  subsequent  chemistry. 

This  work  has  been  published  in  J.  Phys.  Chem.  A,  and  is  available  on  the  web  (DOI: 
10.1021/jp3099889):  Albert  Wagner,  Luis  A  Rivera-Rivera,  Damien  Bachellerie,  Jamin  W 
Perry,  and  Donald  L  Thompson,  "A  Classical  Trajectory  Study  of  the  Dissociation  and 
Isomerization  of  C2H5.”  The  abstract  follows. 

Abstract:  Motivated  by  photo-dissociation  experiments  on  the  ethyl  radical,  we  have 
performed  a  classical  trajectory  study  of  the  dissociation  and  isomerization  of  C2H5  out  to 
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100  ps  over  the  energy  range  100  to  150  kcal/mol.  We  used  a  customized  version  of 
AIREBO  semi-empirical  potential  [Stuart  et  al.  J.  Chem.  Phys.  112,  6472  (2000)]  for  C2H5. 
The  trajectory  results  were  fit  to  a  kinetics  model  of  isomerization/dissociation 
competition  and  the  model  in  turn  was  used  to  project  results  to  the  time  scales  of  the 
experiments  (~100  ns),  and  to  test  for  evidence  of  incomplete  mode-mixing  in  phase  space. 
This  analysis  shows  that  there  is  incomplete  mode  mixing  with  the  AIREBO  potential  at  the 
upper  end  of  our  energy  range,  resulting  in  isomer  population  decay  that  is  poly¬ 
exponential  and  in  dissociation  and  isomerization  rates  that  change  with  the  number  of 
isomerizations  that  have  occurred.  At  the  lower  end  of  our  energy  range,  isomer  population 
decay  is  well  described  by  a  single  exponential  with  isomerization  and  dissociation  rates 
that  are  independent  of  the  number  of  isomerization.  No  evidence  was  found  for  an 
extremely  slow,  non-RRKM  dissociation  process  suggested  by  two  experimental  studies. 
The  fraction  of  energy  deposited  in  translational  motion  of  the  products  is  consistent  with 
previous  theoretical  results  of  others  but  noticeably  lower  than  the  experimental  studies. 
Also  in  contrast  to  the  experimental  and  theoretical  studies  of  others,  these  trajectory 
studies  show  non-negligible  isotopic  scrambling  throughout  our  energy  range  due  to  some 
isomerization  before  dissociation  in  CH3CD2.  The  trajectory  results  support  the  assumption 
that  at  these  energies  an  H  or  a  D  atom  is  equally  likely  to  dissociate  from  the  mixed- 
isotope  methyl  end  of  the  molecule.  We  have  also  carried  out  direct  dynamics  simulations 
for  C2H5  using  the  DFT  method  M05  with  the  6-31+G(d,p)  basis  set  for  total  energies  150, 
200,  and  250  kcal/mol.  These  calculations  are  the  preliminary  step  in  exploring  ways  to 
facilitate  direct  dynamics  with  interpolations  to  decrease  the  number  of  electronic 
structure  calculations. 

The  intramolecular  dynamics  and  unimolecular  dissociation  of  isolated  HO2:  We 

have  carried  out  studies  of  intramolecular  vibrational  energy  transfer,  isomerization,  and 
dissociation  of  HO2  using  highly  accurate  PESs  fit  to  high-level  electronic  structure  theory 
results.  All  the  previous  theoretical  studies  of  the  radical  were  done  using  approximate 
PESs,  and  most  focused  primarily  on  it  as  a  collision  complex  in  reactions  of  H  +  O2  and  OH 
+  0.  We  needed  to  develop  a  better  knowledge  of  the  dynamics  of  the  radical  as  a 
foundation  for  interpreting  our  results  for  pressure  effects.  Our  results  show  that  there  are 
mode-specific  effects  in  both  the  intramolecular  energy  transfer  and  unimolecular 
dissociation  because  the  0-0  and  O-H  vibrational  modes  are  weakly  coupled.  These  results, 
based  on  highly  accurate  PESs,  are  in  conflict  with  conclusions  of  the  previous  theoretical 
studies  that  had  been  carried  out  using  an  approximate  PES.  They  also  illustrate  the  need  to 
develop  more  effective  ways  of  using  high-level  electronic  structure  theory  in  predicting 
vibrational  relaxation  and  chemical  reactions. 

The  intramolecular  dynamics  and  unimolecular  dissociation  of  isolated  HO2:  We 
have  carried  out  studies  of  intramolecular  vibrational  energy  transfer  (IVR),  isomerization, 
and  dissociation  of  HO2  using  highly  accurate  PESs  fit  to  high-level  electronic  structure 
theory  results.  All  the  previous  theoretical  studies  of  the  radical  were  done  using 
approximate  PESs,  and  most  focused  primarily  on  it  as  a  collision  complex  in  reactions  of  H 
+  O2  and  OH  +  0.  We  needed  to  develop  a  better  knowledge  of  the  dynamics  of  the  radical 
as  a  foundation  for  interpreting  our  results  for  pressure  effects.  Our  results  show  that  there 
are  mode-specific  effects  in  both  the  intramolecular  energy  transfer  and  unimolecular 
dissociation  because  the  0-0  and  0-H  vibrational  modes  are  weakly  coupled.  These  results, 
based  on  highly  accurate  PESs,  are  in  conflict  with  conclusions  of  the  previous  theoretical 
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studies  that  had  been  carried  out  using  an  approximate  PES.  The  results  also  illustrate  the 
need  to  develop  more  effective  ways  of  using  high-level  electronic  structure  theory  in 
predicting  vibrational  relaxation  and  chemical  reactions. 

We  have  published  a  study  of  molecular  dynamics  simulations  of  this  important 
combustion  radical  using  highly  accurate  PESs  based  on  high-level  electronic  structure 
theory  results.  The  study  has  been  published:  Jamin  W.  Perry,  Richard  Dawes,  Albert  F. 
Wagner,  and  Donald  L.  Thompson,  "A  Classical  Trajectory  Study  of  the  Intramolecular 
Dynamics,  Isomerization,  and  Unimolecular  Dissociation  of  HO2,"  J.  Chem.  Phys.  139, 
084319  (2013];  http://dx.doi.Org/10.1063/l.4818879j:  the  abstract  follows: 

Abstract:  The  classical  dynamics  and  rates  of  isomerization  and  dissociation  of  HO2 
have  been  studied  using  two  potential  energy  surfaces  (PESs)  based  on  interpolative 
fittings  of  ab  initio  data:  An  interpolative  moving  least-squares  (IMLS)  surface  [A.  Li,  D.  Xie, 
R.  Dawes,  A.  W.  Jasper,  J.  Ma,  and  H.  Guo,  J.  Chem.  Phys.  133 , 144306  (2010)]  and  the  cubic- 
spline-fitted  PES  reported  by  Xu,  Xie,  Zhang,  Lin,  and  Guo  (XXZLG)  [J.  Chem.  Phys.  127, 
024304  (2007)].  Both  PESs  are  based  on  similar,  though  not  identical,  internally  contracted 
multi-reference  configuration  interaction  with  Davidson  correction  (icMRCI+Q)  electronic 
structure  calculations;  the  IMLS  PES  includes  complete  basis  set  (CBS)  extrapolation.  The 
coordinate  range  of  the  IMLS  PES  is  limited  to  non-reactive  processes.  Surfaces-of-section 
show  similar  generally  regular  phase  space  structures  for  the  IMLS  and  XXZLG  PESs  with 
increasing  energy.  The  intramolecular  vibrational  energy  redistribution  (IVR)  at  energies 
above  and  below  the  threshold  of  isomerization  is  slow,  especially  for  0-0  stretch 
excitations,  consistent  with  the  regularity  in  the  surfaces-of-section.  The  slow  IVR  rates 
lead  to  mode-specific  effects  that  are  prominent  for  isomerization  (on  both  the  IMLS  and 
XXZLG)  and  modest  for  unimolecular  dissociation  to  H  +  O2  (accessible  only  on  the  XXZLG 
PES).  Even  with  statistical  distributions  of  initial  energy,  slow  IVR  rates  result  in  double 
exponential  decay  for  isomerization,  with  the  slower  rate  correlated  with  slow  IVR  rates 
for  0-0  vibrational  excitation.  The  IVR  and  isomerization  rates  computed  for  the  IMLS  and 
XXZLG  PESs  are  quantitatively,  but  not  qualitatively,  different  from  one  another  with  the 
largest  differences  ascribed  to  the  ~  2  kcal/mol  difference  in  the  isomerization  barrier 
heights.  The  IMLS  and  XXZLG  results  are  compared  with  those  obtained  using  the  global, 
semi-empirical  double  many-body  expansion  DMBE-IV  PES  [M.  R.  Pastrana,  L.  A.  M. 
Quintales,  J.  Brandao,  and  A.  J.  C.  Varandas,  J.  Chem.  Phys.  94,  8073  (1990)],  for  which  the 
surfaces-of-section  display  more  irregular  phase  space  structure,  much  faster  IVR  rates, 
and  significantly  less  mode-specific  effects  in  isomerization  and  unimolecular  dissociation. 
The  calculated  IVR  results  for  all  three  PESs  are  reasonably  well  represented  by  an  analytic, 
coupled  three-mode  energy  transfer  model. 

Current  Status  on  Improving  ab  initio  PES  fitting  and  direct  dynamics.  In  earlier  work, 
we  developed  reliable,  highly  accurate  methods  for  fitting  high-level  electronic  structure  theory 
results  for  up  to  six  dimensions  (four  atoms).  In  this  grant  we  have  explored  ways  to  speed  up  the 
interpolations  to  obtain  gradients  sufficiently  fast  to  make  IMLS  more  practical  for  larger 
systems.  The  cost  of  evaluations  (interpolations)  for  high-order  basis  functions  can  be 
prohibitively  expensive  in  cpu  time,  thus  limiting  their  use  as  presently  formulated  in  classical 
trajectory  calculations.  We  considered  ways  to  do  this  for  IMLS-directed  global  PES 
construction,  in  trajectory-generated  IMLS  PESs,  and  in  IMLS-assisted  direct  dynamics.  We 
note  however  that  the  IMLS  method,  as  now  formulated,  still  requires  far  too  much  cpu  time  for 
evaluations  of  energies  (and  gradients)  for  effective  use  MD  simulations  of  large  systems.  We 
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have  shown  (in  earlier  studies)  IMLS  to  be  a  powerful  method  for  use  in  predictions  of 
vibrational  spectroscopy  where  far  fewer  energy  points  are  needed.  We  remain  interested  in 
exploring  ways  to  make  better  use  of  high-level  ab  initio  results,  and  thus  will  continue  the  work 
begun  in  this  grant.  Clearly  innovation  is  needed  in  interpolation  methods,  and  other  approaches 
need  to  be  considered.  Another,  perhaps  more  promising,  approach  will  be  improved  functionals 
for  DFT.  We  have  also  carried  out  direct  dynamics  simulations  for  C2H5  using  the  DFT  method 
M05  with  the  6-31+G(d,p)  basis  set  for  total  energies  150,  200,  and  250  kcal/mol. 

Atomic-level  simulations  of  shock  tubes:  We  have  developed  C++  codes  for  the  direct 
atomic-level  simulation  of  a  shock  tube.  The  purpose  is  to  directly  simulate  the  evolution  of 
the  lead  shock  and  secondary  waves  in  the  shock  tube,  including  the  effects  of  reflections  at 
the  ends  of  the  tube  ( e.g .,  re-shocking  of  the  driven  gas,  which  is  what  leads  to  the  high 
temperatures  in  shock  tubes),  and  to  simulate  the  excitation  and  subsequent  relaxation  of 
the  N2  molecules.  In  our  initial  simulations,  the  shock  tube  contained  40,000  He  atoms  as 
the  driver  and  a  10%  mole  ratio  solution  of  N2  in  Ar  as  the  driven  gas  of  55,000  atoms.  The 
initial  temperature  of  the  system  was  298  K.  The  overall  simulation  cell  was  non-periodic 
along  the  length  of  the  shock  tube  but  periodic  in  the  transverse  directions;  5070  A  x  198  A 
x  198  A  ( i.e .,  ~0.5  micron  length  along  the  shock  tube).  The  excitation  and  relaxation  of  N2 
are  being  studied  in  detail  in  terms  of  fundamental  physical  properties  including 
translational  energy,  molecular  angular  momentum,  and  internal  energy.  The  analysis  is 
non-trivial  because  the  system  is  non-steady  in  any  reference  frame;  the  local 
thermodynamic  conditions  under  which  the  N2  excitation  and  relaxation  occurs  vary  along 
the  length  of  the  shock  tube.  Calculations  and  analyses  are  continuing.  We  anticipate 
writing  a  manuscript  for  publication  in  the  next  few  months. 
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